// 版权所有2009 Go作者。保留所有权利。
// 此源代码的使用受BSD样式
// 许可证的约束，该许可证可以在许可证文件中找到。

package math

// 原始C代码和下面的长注释是
// 来自FreeBSD的/usr/src/lib/msun/src/e_sqrt。c和
// 随此通知一同发出。go代码是原始C.
// 
// ==============================================================================================================================
// Sun Microsystems，Inc.版权所有。
// 
// 由Sun Microsystems，Inc.业务部门SunPro开发。
// 使用、复制、修改和分发本
// 软件的权限是免费授予的，前提是保留本通知
// 的权利。
// /===============================================================================================================================
// /
// /\uIEEE754_sqrt（x）
// 返回正确四舍五入的sqrt。
// --------------------------------------
// |如果有硬件，请使用硬件sqrt | 
// ---------------------------------------
// 方法：
// 使用整数算术的逐位方法。（速度慢，但携带方便）
// 1。标准化
// 在[1,4]中用2的偶数幂将x标到y：
// 找到一个整数k，使1<=（y=x*2**（2k））<4，然后
// sqrt（x）=2**k*sqrt（y）
// i0 
// i+12 
// s=2*q，y=2*（y-q）。（1） 
// i 
// 
// 从q计算q，如果（2）为假，则检查
// i+1 i 
// （i+1）
// （i+1）2 
// （q+2）<=y（2）
// i 
// （i+1）
// ，那么q=q；否则q=q+2。
// i+1i+1i 
// 
// ，不难看出，（2）相当于（1+1）
// s+2<=y（3）
// /i 
// 
// （3）的优点是s和y可以通过
// i i 
// 以下递推公式：
// 如果（3）为假
// 
// s=s，y=y；（4） 
// i+1 i+1 i 
// 
// 否则，
// /-i-（i+1）
// /s=s+2，y=y-s-2（5）
// i+1i+1i 
// 
// 人们可以很容易地用归纳法证明（4）和（5）。
// 注意。由于（3）的左侧仅包含i+2位，因此没有必要在（3）中进行完整（53位）比较。
// 3。最终舍入
// 生成53位结果后，我们再计算一位。
// 连同余数，我们可以决定
// 结果是精确的，大于1/2ulp，还是小于1/2ulp 
// （它永远不会等于1/2ulp）。
// 通过检查
// 大+小是否等于大，以及对于某些浮点数“大”和“小”，大-小是否等于大，可以检测舍入模式。
// 
// 
// 注意：舍入模式检测省略。常量“mask”、“shift”、
// 和“bias”在src/math/bits中找到。go 

// Sqrt返回x的平方根。
// 
// 特殊情况是：
// Sqrt（+Inf）=+Inf 
// Sqrt（±0）=±0 
// Sqrt（x<0）=NaN 
// Sqrt（NaN）=NaN 
func Sqrt(x float64) float64 {
	if haveArchSqrt {
		return archSqrt(x)
	}
	return sqrt(x)
}

// 注意：在某些系统上，Sqrt是在汇编中实现的。
// 其他程序集存根跳转到下面的func sqrt。
// 在Sqrt是一条指令的系统上，编译器
// 可以将直接调用转换为直接使用该指令。

func sqrt(x float64) float64 {
	// 特殊情况
	switch {
	case x == 0 || IsNaN(x) || IsInf(x, 1):
		return x
	case x < 0:
		return NaN()
	}
	ix := Float64bits(x)
	// 标准化x 
	exp := int((ix >> shift) & mask)
	if exp == 0 { // 次标准化x 
		for ix&(1<<shift) == 0 {
			ix <<= 1
			exp--
		}
		exp++
	}
	exp -= bias // 无偏指数
	ix &^= mask << shift
	ix |= 1 << shift
	if exp&1 == 1 { // 奇数exp，双x使其偶数
		ix <<= 1
	}
	exp >>= 1 // exp=exp/2，平方根指数
	// 逐位生成sqrt（x）
	ix <<= 1
	var q, s uint64               // q=sqrt（x）
	r := uint64(1 << (shift + 1)) // r=将位从MSB移到LSB 
	for r != 0 {
		t := s + r
		if t <= ix {
			s = t + r
			ix -= t
			q += r
		}
		ix <<= 1
		r >>= 1
	}
	// 最终四舍五入
	if ix != 0 { // 余数，结果不精确
		q += q & 1 // 根据额外位
	}
	ix = q>>1 + uint64(exp-1+bias)<<shift // 有效位+偏置指数
	return Float64frombits(ix)
}
